/* CONSIDER CO-POLLUTANTS FOR BENEFITS */
/* Run lambdas_v1_localpollution.do to generate inputs */ 

clear all
cd "${PATH}/data"


cap log close 
set logtype text
log using "../results/log_PVoutput_analysis3_10KW_lpollution.txt", replace


timer clear
timer on 1
set more off
use "data_15min.dta", clear



// MERGE WITH FILE WITH LAMBDAS (AVOIDED OPERATING COSTS, dAS/dR, AVOIDED EMISSIONS) 
** From file lambdas_v1_localpollution 
preserve
use lambda_15min_locpollution.dta, clear
keep TSO time_utc marginal_cost marginal_emissions marginal_emissions_SOX marginal_emissions_NOX kload solar step_* cost_* emission_* dAS_dR
replace time_utc = time_utc - msofhours(1)
tempfile temp1
save `temp1'
restore

merge m:1 TSO time_utc using `temp1'
keep if _merge == 3 
drop _merge 


// COMPUTE BASELINE OUTPUT
gen year = yofd(dofc(time_utc))


*******************************
// DEFINE LOCALS // 

* CHOOSE incremental step for loop, i.e. XX  MW steps reallocation at the time
local s = $step_size 

* CHOOSE GAMMA // Max share of houses to be covered with solar 
local gamma_min = $gamma_min // 0.05 //
local gamma_max = $gamma_max // 0.50
local gamma_delta = $gamma_delta // 0.05
local number_gammas = round((`gamma_max' - `gamma_min') / `gamma_delta') + 1

* CHOOSE SCC and local pollutants VALUE 
local SCC  31.71  // 100 // 50
local SOX_cost 34051.2 
local NOX_cost 2021.79
*********************************

gen marginal_emissions_EURton = `SCC' * marginal_emissions
gen marginal_emissions_SOX_EURton = `SOX_cost' * marginal_emissions_SOX
gen marginal_emissions_NOX_EURton = `NOX_cost' * marginal_emissions_NOX

// Marginal costs and marginal emissions (mc, me) for each of the steps:
gen mc_step_0 = 0 // Renewables (zero MC)
gen mc_step_1 = 5.1 // Nuclear 
gen mc_step_2 = 6 // Waste and Biomass 
gen mc_step_3 = 11.3  // Lignite
gen mc_step_4 = cost_MWH_coal 
gen mc_step_5 = cost_MWH_ng 
gen mc_step_6 = cost_MWH_crude_oil

gen me_step_0 = 0 // Renewables (zero MC + zero emission)
gen me_step_1 = 0 // Nuclear (no emissions)
gen me_step_2 = 0 // Waste and Biomass (no emissions)
gen me_step_3 = 1.148*`SCC' // Lignite
gen me_step_4 = emission_MWH_coal*`SCC' // hard coal
gen me_step_5 = emission_MWH_ng*`SCC' // gas
gen me_step_6 = emission_MWH_crude_oil*`SCC' // oil

// LOCAL POLLUTANTS [see EU-emissions-data-localpollution.do]
gen me_step_0_locpol = 0 // Renewables (zero MC + zero emission)
gen me_step_1_locpol = 0 // Nuclear (no emissions)
gen me_step_2_locpol = 0 // Waste and Biomass (no emissions)
gen me_step_3_locpol = .0007348*`SOX_cost' + .00085*`NOX_cost'  // Lignite 
gen me_step_4_locpol = .0004455*`SOX_cost' + .0006804*`NOX_cost'  // hard coal
gen me_step_5_locpol = .0002644*`NOX_cost'  // gas
gen me_step_6_locpol = .0003164*`NOX_cost'  // oil

drop cost_* emission_* solar


*Update information on housing stock
gen 	buildings = 3518109 if TSO == "50Hertz"
replace buildings = 6079790 if TSO == "Amprion"
replace buildings = 7034498 if TSO == "TenneT"
replace buildings = 2240361 if TSO == "TransnetBW"

local avg_inst_size = 6.7
gen max_cap = (buildings*`avg_inst_size') /1000 // max cap in MW if all houses had avg. sized solar installation
gen max_cap_gamma = . // max_cap * gamma (in loop)

keep if solar_gen >0 & !missing(solar_gen) // LIMIT sample only to positive solar output hours
****************

** Fraction of small term solar in total production ** 

** 2016 (use max data)
** (FROM solar_capacityTSO.dta, < 10 KW, total cum):
local solarcap_10kw_50Hertz   536.55963
local solarcap_10kw_Amprion   1779.9152
local solarcap_10kw_TenneT    2293.709
local solarcap_10kw_TransnetBW  1117.9264

local fraction_50Hertz   `solarcap_10kw_50Hertz' /   9842.2686
local fraction_Amprion   `solarcap_10kw_Amprion' /   9549.0361
local fraction_TenneT  	 `solarcap_10kw_TenneT'  /   15509.563
local fraction_TransnetBW `solarcap_10kw_TransnetBW' /   5605.6885

// S IS THE SUM OF THE solarcap_10kw_*
local S = 5728 // 5728.1102  // max values of installed capacity (<10KW) in 2016

// Production of residential installations in MWh
rename solar_gen solar_gen_official
gen solar_gen  = .
gen solar_gen_MW = .

levelsof TSO, local(T)
foreach t in `T' {
	replace solar_gen = solar_gen_official * `fraction_`t'' if TSO =="`t'" // production of residential solar systems 
	replace solar_gen_MW = solar_gen_official *`fraction_`t'' /  `solarcap_10kw_`t'' if TSO =="`t'" // production of res. systems [per MW]
}

// For reallocation: solar generation per capacity that is moved around // 
gen solar_gen_s = solar_gen_MW * `s'

*********


// S.D.s FROM LOAD AND SOLAR_SCALED (FROM load_regressions.do)

local factor_sd = 2

local sd_load_50Hertz `factor_sd' * 327.2129
local sd_load_Amprion `factor_sd' * 193.0014
local sd_load_TenneT `factor_sd' * 159.003
local sd_load_TransnetBW `factor_sd' * 93.40471

** Original
local sd_solar_50Hertz `factor_sd' * .0000157
local sd_solar_Amprion `factor_sd' * .0000192
local sd_solar_TenneT `factor_sd' * .000028
local sd_solar_TransnetBW `factor_sd' * .0000219

gen MB_bench  = -dAS_dR + marginal_cost +  marginal_emissions_EURton + marginal_emissions_SOX_EURton + marginal_emissions_NOX_EURton

*Calculate total value: MB*solar output
gen MB_bench_tot = MB_bench * solar_gen
gen MB_bench_AS = -dAS_dR * solar_gen
gen MB_bench_MC = marginal_cost * solar_gen
gen MB_bench_ME = marginal_emissions_EURton * solar_gen
gen MB_bench_ME_locpol = (marginal_emissions_SOX_EURton + marginal_emissions_NOX_EURton) * solar_gen


// AGGREGATE TOTAL BENEFITS 
preserve

	// Collapse by TSO
	collapse (sum) MB_bench_* , by(TSO)
	
	qui: su MB_bench_tot  if TSO == "50Hertz"
		local value_50Hertz_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "Amprion"	
		local value_Amprion_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TenneT"	
		local value_TenneT_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TransnetBW"
		local value_TransnetBW_bench = r(mean)
 
	qui: su MB_bench_AS  if TSO == "50Hertz"
		local MB_bench_AS_50Hertz = r(mean)
	qui: su MB_bench_AS  if TSO == "Amprion"
		local MB_bench_AS_Amprion = r(mean)
	qui: su MB_bench_AS  if TSO == "TenneT"
		local MB_bench_AS_TenneT = r(mean)
	qui: su MB_bench_AS  if TSO == "TransnetBW"
		local MB_bench_AS_TransnetBW = r(mean)

	qui: su MB_bench_MC  if TSO == "50Hertz"
		local MB_bench_MC_50Hertz = r(mean)
	qui: su MB_bench_MC  if TSO == "Amprion"
		local MB_bench_MC_Amprion = r(mean)
	qui: su MB_bench_MC  if TSO == "TenneT"
		local MB_bench_MC_TenneT = r(mean)
	qui: su MB_bench_MC  if TSO == "TransnetBW"
		local MB_bench_MC_TransnetBW = r(mean)
		
	qui: su MB_bench_ME  if TSO == "50Hertz"
		local MB_bench_ME_50Hertz = r(mean)
	qui: su MB_bench_ME  if TSO == "Amprion"
		local MB_bench_ME_Amprion = r(mean)
	qui: su MB_bench_ME  if TSO == "TenneT"
		local MB_bench_ME_TenneT = r(mean)
	qui: su MB_bench_ME  if TSO == "TransnetBW"
		local MB_bench_ME_TransnetBW = r(mean)

	qui: su MB_bench_ME_locpol  if TSO == "50Hertz"
		local MB_bench_ME_locpol_50Hertz = r(mean)
	qui: su MB_bench_ME_locpol  if TSO == "Amprion"
		local MB_bench_ME_locpol_Amprion = r(mean)
	qui: su MB_bench_ME_locpol  if TSO == "TenneT"
		local MB_bench_ME_locpol_TenneT = r(mean)
	qui: su MB_bench_ME_locpol  if TSO == "TransnetBW"
		local MB_bench_ME_locpol_TransnetBW = r(mean)
	
	
	// Collapse for total value
	collapse (sum) MB_bench_* 
	qui: su MB_bench_tot
		local value_bench_total = r(mean)
	qui: su MB_bench_AS
		local value_bench_AS = r(mean)
	qui: su MB_bench_MC
		local value_bench_MC = r(mean)
	qui: su MB_bench_ME
		local value_bench_ME = r(mean)
	qui: su MB_bench_ME_locpol
		local value_bench_ME_locpol = r(mean)
restore

drop MB_bench MB_bench_* dAS_dR 

*********
// Generate Matrices to store values
*Matrix 1: Allocation, Value. by TSO, total
mat ratios = J(`number_gammas', 20, .)

*Matrix 2: Decomposition of gains: ASC, MC, ME. by TSO, total
mat ratios_decomposition = J(`number_gammas', 42, .)

mat ratios_sd_plus_1_plus_1 = J(`number_gammas', 20, .)
mat ratios_sd_minus_1_minus_1 = J(`number_gammas', 20, .)
mat ratios_sd_minus_1_plus_1 = J(`number_gammas', 20, .)
mat ratios_sd_plus_1_minus_1 = J(`number_gammas', 20, .)

*********
// Reallocation //

cap erase temp_v4.dta
save temp_v4.dta

forvalues bands = 0 / 4 {
** DEBUG
//local bands = 0 
	
	use temp_v4, clear
	
	if `bands' == 1 {
		// CASE + +
		local sign_load  1
		local sign_solar 1
	}
	else if `bands' == 2 {
		// CASE - - 
		local sign_load  -1
		local sign_solar -1
	}
	else if `bands' == 3 {
		// CASE + -
		local sign_load  1 
		local sign_solar -1		
	}
	else if `bands' == 4 {
		// CASE - +
		local sign_load  -1
		local sign_solar 1 		
	}

	
	if `bands' > 0 {
		replace load = load + `sign_load' * `sd_load_50Hertz' if TSO == "50Hertz"
		replace load = load + `sign_load' * `sd_load_Amprion' if TSO == "Amprion"
		replace load = load + `sign_load' * `sd_load_TenneT' if TSO == "TenneT"
		replace load = load + `sign_load' * `sd_load_TransnetBW' if TSO == "TransnetBW"
		
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_50Hertz' > 0, solar_gen_MW + `sign_solar' * `sd_solar_50Hertz', 0) if TSO == "50Hertz"
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_Amprion' > 0, solar_gen_MW + `sign_solar' * `sd_solar_Amprion', 0) if TSO == "Amprion"
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_TenneT' > 0, solar_gen_MW + `sign_solar' * `sd_solar_TenneT', 0) if TSO == "TenneT"
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_TransnetBW' > 0, solar_gen_MW + `sign_solar' * `sd_solar_TransnetBW', 0) if TSO == "TransnetBW"
	
	// NEED TO REDEFINE solar_gen_s
		replace solar_gen_s = solar_gen_MW * `s'
	
	// Check for all load to be positive
		replace load = 0 if load < 0  // in 25 out of 149054 obs. 
	}

	

	local j 1
	* Need to calculate the value for each MW of installed solar and store allocation
	forvalues gamma = `gamma_min'(`gamma_delta')`gamma_max' { 	// LOOP OVER GAMMA VALUES

	
		*Generate variable to keep track of installed in reallocation  
		cap drop solar_cap solar_capl1  cum_solar_gen solar_new value_50Hertz  value_Amprion value_TenneT value_TransnetBW count_excess_load
		cap drop value_AS_* value_MC_* value_ME_*
		gen solar_cap = 0 // solar capacity in the loop (i)
		gen solar_capl1 = 0 // solar capacity at iteration (i-1)
		gen cum_solar_gen = 0 // keep track of cumulative solar production -> for correct "step" on supply curve
		gen solar_new = . // define difference of solar_reallocated - solar_original 
		
		gen count_excess_load = 0 // keep track how many times cumulative solar > load
		
		// Store values: Total and decomposition: ASC, MC, ME
		gen value_50Hertz = 0
		gen value_Amprion = 0
		gen value_TenneT = 0
		gen value_TransnetBW = 0
		
		gen value_AS_50Hertz = 0
		gen value_AS_Amprion = 0
		gen value_AS_TenneT = 0
		gen value_AS_TransnetBW = 0
		gen value_MC_50Hertz = 0
		gen value_MC_Amprion = 0
		gen value_MC_TenneT = 0
		gen value_MC_TransnetBW = 0
		gen value_ME_50Hertz = 0
		gen value_ME_Amprion = 0
		gen value_ME_TenneT = 0
		gen value_ME_TransnetBW = 0
		gen value_ME_locpol_50Hertz = 0
		gen value_ME_locpol_Amprion = 0
		gen value_ME_locpol_TenneT = 0
		gen value_ME_locpol_TransnetBW = 0
		
		// FREQUENCY OF TIMES TO THE RIGHT OR TO THE LEFT OF INITIAL EQM POINT
		local in_flat_MB 0
		local in_step_6 0
		local in_step_5 0
		local in_step_4 0
		local in_step_3 0
		local in_step_2 0
		local in_step_1 0
		local in_step_0 0
		
		
		forvalues i = `s'(`s')`S' {  // loop over all MW of installed capacity to determine optimal allocation
			
			di " " 
			di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% " 
			di "Reallocating MW: `i' of `S'"
		
			// Set max. capacity per TSO 
			replace max_cap_gamma = max_cap * `gamma'
		
		
			// Cumulative solar generation: include additional MW of solar installed
			replace solar_cap = solar_cap + `s'
			replace cum_solar_gen =  solar_gen_MW * solar_cap
			* solar_gen_s calculated in line 160

		// DERIVATIVES OF ABOVE REGRESSION MODEL WRT SOLAR
		gen dAS_dR = 5.174 +  (-0.00223) * 2 * solar_gen + ///
			( -9.83e-09) * 3 * solar_gen^2  + (0.000380) * load + ///
			 ( -6.76e-08) * load^2 + ( 0.000000189) * 2 * solar_gen * load if kload == 1 & TSO == "50Hertz"
			
		replace dAS_dR =  6.455 +  (-0.000391) * 2 * solar_gen + ///
			(-3.10e-08) * 3 * solar_gen^2  + ( 0.0000577) * load + ///
			 (-6.62e-08) * load^2 + (5.46e-08) * 2 * solar_gen * load if kload == 3 & TSO == "50Hertz"	
			
		replace dAS_dR =   18.44 +   (-0.00560) * 2 * solar_gen + ///
			(0.000000512) * 3 * solar_gen^2  + ( -0.00260) * load + ///
			  0.000000159 * load^2 + (0.000000185) * 2 * solar_gen * load if kload == 2 & TSO == "50Hertz"			  

		replace dAS_dR = 79.60 +  0.00101 * 2 * solar_gen + ///
			(0.000000101) * 3 * solar_gen^2  + (-0.00558) * load + ///
			 0.000000103 * load^2 + (-9.49e-08) * 2 * solar_gen * load if kload == 1 & TSO == "Amprion"
			
		replace dAS_dR =  -20.40 +  (-0.000250) * 2 * solar_gen + ///
			(0.000000442) * 3 * solar_gen^2  + (0.00184) * load + ///
			 -2.90e-08 * load^2 + (-0.000000117) * 2 * solar_gen * load if kload == 2 & TSO == "Amprion"	
			
		replace dAS_dR =  77.84 + 0.00388 * 2 * solar_gen + ///
			(-6.73e-09) * 3 * solar_gen^2  + (-0.00705) * load + ///
			  0.000000156 * load^2 + (-0.000000155) * 2 * solar_gen * load if kload == 3 & TSO == "Amprion"			  

		replace dAS_dR = 380.1 + (-0.000817) * 2 * solar_gen + ///
			(1.43e-08) * 3 * solar_gen^2  + (-0.0347) * load + ///
			  0.000000788 * load^2 + (3.51e-08) * 2 * solar_gen * load if kload == 1 & TSO == "TenneT"
			
		replace dAS_dR = -3.852 + 0.00105 * 2 * solar_gen + ///
			(-5.31e-09) * 3 * solar_gen^2  + (-0.0000909) * load + ///
			 1.98e-08 * load^2 + (-5.87e-08) * 2 * solar_gen * load if kload == 2 & TSO == "TenneT"	
			
		replace dAS_dR = 42.15 + 0.00161 * 2 * solar_gen + ///
			(-1.91e-08) * 3 * solar_gen^2  + (-0.00512) * load + ///
			  0.000000148 * load^2 + (-6.88e-08) * 2 * solar_gen * load if kload == 3 & TSO == "TenneT"			  

		replace dAS_dR = 130.8 +  0.00561 * 2 * solar_gen + ///
			(-1.23e-08) * 3 * solar_gen^2  + (-0.0329) * load + ///
			 0.00000202 * load^2 + (-0.000000616) * 2 * solar_gen * load if kload == 1 & TSO == "TransnetBW"
			
		replace dAS_dR = 40.25 + 0.0160 * 2 * solar_gen + ///
			(0.00000132) * 3 * solar_gen^2  + (-0.0217) * load + ///
			 0.00000274 * load^2 + (-0.00000375) * 2 * solar_gen * load if kload == 2 & TSO == "TransnetBW"	
			
		replace dAS_dR = 105.7 + 0.0130 * 2 * solar_gen + ///
			(-0.000000160) * 3 * solar_gen^2  + (-0.0314) * load + ///
			  0.00000223 * load^2 + (-0.00000150) * 2 * solar_gen * load if kload == 3 & TSO == "TransnetBW"			   
				
		
		
		// Define location on supply curve -> identify MC and marginal Emissions
			replace solar_new =  cum_solar_gen - solar_gen // how much more solar do I generate compared to original allocation? 
			
		
			if (cum_solar_gen <= solar_gen) { 
				*Apply E[MB] as "observed": assumption of "flat MB curves"
				gen MB = -dAS_dR + marginal_cost + marginal_emissions_EURton + marginal_emissions_SOX_EURton + marginal_emissions_NOX_EURton
				** Track how many times we are in flat part 
				gen MB_tot = cond(cum_solar_gen < load, MB * solar_gen_s, 0) // evaluate MB for block of reallocated solar
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0)
				gen MB_MC = cond(cum_solar_gen  < load, marginal_cost * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, marginal_emissions_EURton * solar_gen_s , 0)
				gen MB_ME_locpol = cond(cum_solar_gen  < load, (marginal_emissions_SOX_EURton + marginal_emissions_NOX_EURton) * solar_gen_s , 0)
		
				local in_flat_MB `in_flat_MB' + 1
			}
		
		
			else {
				*We reallocate more solar to TSO as in "benchmark" case -> need to adjust benefits
			
				replace count_excess_load = count_excess_load + 1 // One more interation for which cum solar output > load
			
				if ((solar_new > 0) & (solar_new <= step_6)) {
					* Apply MB from step6
					gen MB = -dAS_dR + mc_step_6 + me_step_6 + me_step_6_locpol
					
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0)
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_6 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_6 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_6_locpol * solar_gen_s, 0)
					
					local in_step_6 `in_step_6' + 1
		
				}
		
				else if ((solar_new > step_6) & (solar_new <= (step_5 + step_6))) {
					*Apply MB from step 5
					gen MB = -dAS_dR + mc_step_5 + me_step_5 + me_step_5_locpol
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_5 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_5 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_5_locpol * solar_gen_s, 0)
					
					local in_step_5 `in_step_5' + 1
				}
		
				else if  ((solar_new > step_6 + step_5) & (solar_new <= (step_4 + step_5 + step_6))) {
					* Apply MB from step 4
					gen MB = -dAS_dR + mc_step_4 + me_step_4 + me_step_4_locpol
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_4 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_4 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_4_locpol * solar_gen_s, 0)
					
					local in_step_4 `in_step_4' + 1
				}
		
				else if  ((solar_new > step_6 + step_5 + step_4) & (solar_new <= (step_3 + step_4 + step_5 + step_6))) {
					* Apply MB from step 3
					gen MB = -dAS_dR + mc_step_3 + me_step_3 + me_step_3_locpol
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_3 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_3 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_3_locpol * solar_gen_s, 0)

					local in_step_3 `in_step_3' + 1
				}
		
				else if  ((solar_new > step_6 + step_5 + step_4 + step_3) & (solar_new <= (step_2 + step_3 + step_4 + step_5 + step_6))) {
					* Apply MB from step 2
					gen MB = -dAS_dR + mc_step_2 + me_step_2 + me_step_2_locpol
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_2 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_2 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_2_locpol * solar_gen_s, 0)
					
					local in_step_2 `in_step_2' + 1
				}
		
				else if  ((solar_new > step_6 + step_5 + step_4 + step_3 + step_2) & (solar_new <= (step_1 + step_2 + step_3 + step_4 + step_5 + step_6))) {
					* Apply MB from step 1
					gen MB = -dAS_dR + mc_step_1 + me_step_1 + me_step_1_locpol
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0)
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_1 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_1 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_1_locpol * solar_gen_s, 0)

					local in_step_1 `in_step_1' + 1
				}
		
				else {
					* Apply MB from step 0
					gen MB = -dAS_dR + mc_step_0 + me_step_0 + me_step_0_locpol
					gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen  < load, mc_step_0 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen  < load, me_step_0 * solar_gen_s, 0)
					gen MB_ME_locpol = cond(cum_solar_gen  < load, me_step_0_locpol * solar_gen_s, 0)

					local in_step_0 `in_step_0' + 1
				}
		
			} //else
		
		
			** OBTAIN RANKING IN AUTOMATED FASHION ** 
			egen meanMB = mean(-MB_tot), by(TSO)
			egen MB_order = axis(meanMB TSO), label(TSO)
			qui: tabstat MB_tot , by(MB_order) s(mean sd N) save
			local first   `r(name1)' 
			local second  `r(name2)' 
			local third   `r(name3)' 
			local fourth  `r(name4)'
			* su MB if TSO == "`first'"
			di "RANKING TSO: `first', `second', `third', `fourth'"
	
		
			// Allocation rule: nest options // 
			** Store total solar installed, cumulative value of solar installed
		
			gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`first'" 
			replace alloc_constraint = 0 if alloc_constraint!=1 // replaces value in all TSOs
			qui: su alloc_constraint
			local alloc_1 =  `r(mean)'
		
			if `alloc_1' == 0 {
				replace solar_cap = solar_cap - `s' if (TSO == "`second'" | TSO == "`third'" | TSO == "`fourth'")
				* SET back to original value for the TSO's where no allocation - consider constaints!! 
			}
		
		
			else if `alloc_1' != 0 {
		
				drop alloc_constraint
				gen alloc_constraint = 1 if  (solar_cap > max_cap_gamma) & TSO == "`second'" 
				replace alloc_constraint = 0 if alloc_constraint!=1
				qui: su alloc_constraint
				local alloc_2 =  `r(mean)'
		
				if `alloc_2' == 0 {
					replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`third'" | TSO == "`fourth'") 
					* SET back to original value for the TSO's where no allocation - consider constaints!! 
				}
		
				else if `alloc_2' != 0 {
		
					drop alloc_constraint
					gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`third'" 
					replace alloc_constraint = 0 if alloc_constraint!=1
					qui: su alloc_constraint
					local alloc_3 =  `r(mean)'
		
					if `alloc_3' == 0 {
						replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`fourth'") 
						* SET back to original value for the TSO's where no allocation - consider constraints!! 
					}
		
					else {
						// Put in forth TSO 
						replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`third'") 
					} // TSO 4
				} // TSO 3
			} // TSO 2
		
			// Store total value and in case "decomposition"
			qui: levelsof TSO, local(T)
			foreach t in `T' {
				replace value_`t' = value_`t' + MB_tot if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_AS_`t' = value_AS_`t' + MB_AS if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_MC_`t' = value_MC_`t' + MB_MC if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_ME_`t' = value_ME_`t' + MB_ME if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_ME_locpol_`t' = value_ME_locpol_`t' + MB_ME_locpol if (solar_cap != solar_capl1) & TSO == "`t'"

				}
			// value_TSO: sums total MB for each block of solar that is reallocated
		
			replace solar_capl1 = solar_cap
			drop dAS_dR alloc_constraint MB_order meanMB MB MB_tot MB_AS MB_MC MB_ME MB_ME_locpol
			
		
		} // loop over S (total capacity to reallocate)
		
		
		// AGGREGATE TOTAL BENEFITS 
		qui: su solar_cap if TSO =="50Hertz"
		local solar_cap_50Hertz = r(mean) 
		
		qui: su solar_cap if TSO =="Amprion"
		local solar_cap_Amprion = r(mean)
		
		qui: su solar_cap if TSO =="TenneT"
		local solar_cap_TenneT = r(mean)
		
		qui: su solar_cap if TSO =="TransnetBW"
		local solar_cap_TransnetBW = r(mean)
		
			
		
		// Value by TSO
		
		preserve
		collapse (sum) value* 
		
		qui: su value_50Hertz
			local value_50Hertz = r(mean)
		qui: su value_Amprion
			local value_Amprion = r(mean)
		qui: su value_TenneT
			local value_TenneT = r(mean)
		qui: su  value_TransnetBW
			local value_TransnetBW = r(mean)
			
		// Total of average MB over two year period
		local value_total = `value_50Hertz' + `value_Amprion' + `value_TenneT' + `value_TransnetBW'		
		
		qui: su value_AS_50Hertz  
			local MB_AS_50Hertz = r(mean)
		qui: su value_AS_Amprion  
			local MB_AS_Amprion = r(mean)
		qui: su value_AS_TenneT
			local MB_AS_TenneT = r(mean)
		qui: su value_AS_TransnetBW
			local MB_AS_TransnetBW = r(mean)
		
		qui: su value_MC_50Hertz  
			local MB_MC_50Hertz = r(mean)
		qui: su value_MC_Amprion  
			local MB_MC_Amprion = r(mean)
		qui: su value_MC_TenneT
			local MB_MC_TenneT = r(mean)
		qui: su value_MC_TransnetBW
			local MB_MC_TransnetBW = r(mean)		
		
		qui: su value_ME_50Hertz  
			local MB_ME_50Hertz = r(mean)
		qui: su value_ME_Amprion  
			local MB_ME_Amprion = r(mean)
		qui: su value_ME_TenneT
			local MB_ME_TenneT = r(mean)
		qui: su value_ME_TransnetBW
			local MB_ME_TransnetBW	= r(mean)	

		qui: su value_ME_locpol_50Hertz  
			local MB_ME_locpol_50Hertz = r(mean)
		qui: su value_ME_locpol_Amprion  
			local MB_ME_locpol_Amprion = r(mean)
		qui: su value_ME_locpol_TenneT
			local MB_ME_locpol_TenneT = r(mean)
		qui: su value_ME_locpol_TransnetBW
			local MB_ME_locpol_TransnetBW	= r(mean)	
		
		egen value_AS = rowtotal(value_AS_*)
		egen value_MC = rowtotal(value_MC_*)
		egen value_ME = rowtotal(value_ME_50Hertz value_ME_Amprion value_ME_TenneT value_ME_TransnetBW) // need to adjust here, otherwise double counting of emissions benefits
		egen value_ME_locpol = rowtotal(value_ME_locpol_*)
		
		qui: su value_AS
			local value_AS_tot = r(mean)
		
		qui: su value_MC
			local value_MC_tot = r(mean)
		
		qui: su value_ME
			local value_ME_tot = r(mean)	

		qui: su value_ME_locpol
			local value_ME_locpol_tot = r(mean)	
						
		restore	
			
		di "Bands is `bands'"
		di " "
		
		// Fill in Matrices
	
		if `bands' == 0 {
			* Mat 1: Gains
			mat ratios[`j',1] = `SCC'
			mat ratios[`j',2] = `gamma'
			mat ratios[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios[`j',4] = `solar_cap_50Hertz' 
			mat ratios[`j',5] = `solarcap_10kw_Amprion'
			mat ratios[`j',6] = `solar_cap_Amprion'
			mat ratios[`j',7] = `solarcap_10kw_TenneT'
			mat ratios[`j',8] = `solar_cap_TenneT'
			mat ratios[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios[`j',10] = `solar_cap_TransnetBW' 
			mat ratios[`j',11] = `value_bench_total' / 1e6
			mat ratios[`j',12] = `value_total' / 1e6
			mat ratios[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios[`j',14] = `value_50Hertz' / 1e6
			mat ratios[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios[`j',16] = `value_Amprion' / 1e6
			mat ratios[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios[`j',18] = `value_TenneT' / 1e6
			mat ratios[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios[`j',20] = `value_TransnetBW' / 1e6
				
			*Mat 2: Decomposition
			mat ratios_decomposition[`j',1] = `SCC'
			mat ratios_decomposition[`j',2] = `gamma'
			mat ratios_decomposition[`j',3] = `value_bench_AS' / 1e6
			mat ratios_decomposition[`j',4] = `value_AS_tot' / 1e6
			mat ratios_decomposition[`j',5] = `value_bench_MC' / 1e6
			mat ratios_decomposition[`j',6] =  `value_MC_tot' / 1e6
			mat ratios_decomposition[`j',7] = `value_bench_ME' / 1e6
			mat ratios_decomposition[`j',8] = `value_ME_tot' / 1e6
			mat ratios_decomposition[`j',9] = `value_bench_ME_locpol' / 1e6
			mat ratios_decomposition[`j',10] = `value_ME_locpol_tot' / 1e6		
			mat ratios_decomposition[`j',11] = `MB_bench_AS_50Hertz' / 1e6
			mat ratios_decomposition[`j',12] = `MB_AS_50Hertz' / 1e6
			mat ratios_decomposition[`j',13] = `MB_bench_MC_50Hertz' / 1e6
			mat ratios_decomposition[`j',14] = `MB_MC_50Hertz' / 1e6
			mat ratios_decomposition[`j',15] = `MB_bench_ME_50Hertz' / 1e6
			mat ratios_decomposition[`j',16] = `MB_ME_50Hertz' / 1e6
			mat ratios_decomposition[`j',17] = `MB_bench_ME_locpol_50Hertz' / 1e6
			mat ratios_decomposition[`j',18] = `MB_ME_locpol_50Hertz' / 1e6
			mat ratios_decomposition[`j',19] = `MB_bench_AS_Amprion' / 1e6
			mat ratios_decomposition[`j',20] = `MB_AS_Amprion' / 1e6
			mat ratios_decomposition[`j',21] = `MB_bench_MC_Amprion' / 1e6
			mat ratios_decomposition[`j',22] = `MB_MC_Amprion' / 1e6
			mat ratios_decomposition[`j',23] = `MB_bench_ME_Amprion' / 1e6
			mat ratios_decomposition[`j',24] = `MB_ME_Amprion' / 1e6
			mat ratios_decomposition[`j',25] = `MB_bench_ME_locpol_Amprion' / 1e6
			mat ratios_decomposition[`j',26] = `MB_ME_locpol_Amprion' / 1e6
			mat ratios_decomposition[`j',27] = `MB_bench_AS_TenneT' / 1e6
			mat ratios_decomposition[`j',28] = `MB_AS_TenneT' / 1e6
			mat ratios_decomposition[`j',29] = `MB_bench_MC_TenneT' / 1e6
			mat ratios_decomposition[`j',30] = `MB_MC_TenneT' / 1e6
			mat ratios_decomposition[`j',31] = `MB_bench_ME_TenneT' / 1e6
			mat ratios_decomposition[`j',32] = `MB_ME_TenneT' / 1e6
			mat ratios_decomposition[`j',33] = `MB_bench_ME_locpol_TenneT' / 1e6
			mat ratios_decomposition[`j',34] = `MB_ME_locpol_TenneT' / 1e6
			mat ratios_decomposition[`j',35] = `MB_bench_AS_TransnetBW' / 1e6
			mat ratios_decomposition[`j',36] = `MB_AS_TransnetBW' / 1e6
			mat ratios_decomposition[`j',37] = `MB_bench_MC_TransnetBW' / 1e6
			mat ratios_decomposition[`j',38] = `MB_MC_TransnetBW' / 1e6
			mat ratios_decomposition[`j',39] = `MB_bench_ME_TransnetBW' / 1e6
			mat ratios_decomposition[`j',40] = `MB_ME_TransnetBW' / 1e6
			mat ratios_decomposition[`j',41] = `MB_bench_ME_locpol_TransnetBW' / 1e6
			mat ratios_decomposition[`j',42] = `MB_ME_locpol_TransnetBW' / 1e6			
			** ALL VALUES (value, MB, expressed in Million EUR)
		}	
		else if `bands' == 1 {
			* Mat 3: Gains 
			mat ratios_sd_plus_1_plus_1[`j',1] = `SCC'
			mat ratios_sd_plus_1_plus_1[`j',2] = `gamma'
			mat ratios_sd_plus_1_plus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_plus_1_plus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_plus_1_plus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_plus_1_plus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_plus_1_plus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_plus_1_plus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_plus_1_plus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_plus_1_plus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_plus_1_plus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',20] = `value_TransnetBW' / 1e6
		}
		else if `bands' == 2 {
			* Mat 4: Gains 
			mat ratios_sd_minus_1_minus_1[`j',1] = `SCC'
			mat ratios_sd_minus_1_minus_1[`j',2] = `gamma'
			mat ratios_sd_minus_1_minus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_minus_1_minus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_minus_1_minus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_minus_1_minus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_minus_1_minus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_minus_1_minus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_minus_1_minus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_minus_1_minus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_minus_1_minus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',20] = `value_TransnetBW' / 1e6
		}
		else if `bands' == 3 {
			* Mat 4: Gains 
			mat ratios_sd_plus_1_minus_1[`j',1] = `SCC'
			mat ratios_sd_plus_1_minus_1[`j',2] = `gamma'
			mat ratios_sd_plus_1_minus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_plus_1_minus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_plus_1_minus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_plus_1_minus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_plus_1_minus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_plus_1_minus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_plus_1_minus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_plus_1_minus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_plus_1_minus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',20] = `value_TransnetBW' / 1e6
		}
		else if `bands' == 4 {
			* Mat 4: Gains 
			mat ratios_sd_minus_1_plus_1[`j',1] = `SCC'
			mat ratios_sd_minus_1_plus_1[`j',2] = `gamma'
			mat ratios_sd_minus_1_plus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_minus_1_plus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_minus_1_plus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_minus_1_plus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_minus_1_plus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_minus_1_plus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_minus_1_plus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_minus_1_plus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_minus_1_plus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',20] = `value_TransnetBW' / 1e6
		}
					
		// CHECK LOCATION ON SUPPLY CURVE
		di " "
		di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "
		di " "
		di "FREQUENCY OF TIMES TO THE LEFT AND TO THE RIGHT OF INITIAL EQM POINT"
		di `in_flat_MB'
		di `in_step_6'
		di `in_step_5'
		di `in_step_4'
		di `in_step_3'
		di `in_step_2'
		di `in_step_1'
		di `in_step_0'
			
		di " "
		di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "
		di "Number of times solar production is larger than load "
		qui: su count_excess_load if TSO == "50Hertz"
		di "50 Hertz: " `r(mean)'
		qui: su count_excess_load if TSO == "TenneT"
		di "TenneT: " `r(mean)'
		qui: su count_excess_load if TSO == "Amprion"
		di "Amprion: " `r(mean)'
		qui: su count_excess_load if TSO == "TransnetBW"
		di "TransnetBW: " `r(mean)'
		
		
		local ++j
		di `j'

	} // LOOP OVER GAMMA

} // LOOP OVER BANDS	



// Extract values + plot
mat colnames ratios = SCC gamma cap_50HZ cap_50HZ_alloc cap_Amprion cap_Amprion_alloc cap_TenneT cap_TenneT_alloc cap_Transnet cap_Transnet_alloc v_bench v_alloc /// 
v_50HZ v_50HZ_alloc v_Amprion v_Amprion_alloc v_TenneT v_TenneT_alloc v_Transnet v_Transnet_alloc
svmat ratios, names(col)
	
mat colnames ratios_decomposition = SCC_ gamma_ v_AS v_AS_alloc v_MC v_MC_alloc v_ME v_ME_alloc v_ME_locpol v_ME_locpol_alloc v_AS_50HZ v_AS_50HZ_alloc v_MC_50HZ v_MC_50HZ_alloc ///
v_ME_50HZ v_ME_50HZ_alloc v_ME_locpol_50HZ v_ME_locopol_50HZ_alloc v_AS_Amp v_AS_50_HZ_alloc v_MC_Amp v_MC_Amp_alloc v_ME_Amp v_ME_Amp_alloc ///
v_ME_locpol_Amp v_ME_locpol_Amp_alloc  v_AS_TenneT v_AS_TenneT_alloc v_MC_TenneT v_MC_TenneT_alloc v_ME_TenneT v_ME_TenneT_alloc v_ME_locpol_TenneT v_ME_locpol_TenneT_alloc v_AS_Trans v_AS_Trans_alloc v_MC_Trans v_MC_Trans_alloc v_ME_Trans v_ME_Trans_alloc v_ME_locpol_Trans v_ME_locpol_Trans_alloc
svmat ratios_decomposition, names(col)

		
mat colnames ratios_sd_plus_1_plus_1 = SCC_sd_p_1_p_1 gamma_sd_p_1_p_1 cap_50HZ_sd_p_1_p_1 cap_50HZ_alloc_sd_p_1_p_1 cap_Amprion_sd_p_1_p_1 cap_Amprion_alloc_sd_p_1_p_1 cap_TenneT_sd_p_1_p_1 cap_TenneT_alloc_sd_p_1_p_1 cap_Transnet_sd_p_1_p_1 cap_Transnet_alloc_sd_p_1_p_1 v_bench_sd_p_1_p_1 v_alloc_sd_p_1_p_1 /// 
v_50HZ_sd_p_1_p_1 v_50HZ_alloc_sd_p_1_p_1 v_Amprion_sd_p_1_p_1 v_Amprion_alloc_sd_p_1_p_1 v_TenneT_sd_p_1_p_1 v_TenneT_alloc_sd_p_1_p_1 v_Transnet_sd_p_1_p_1 v_Transnet_alloc_sd_p_1_p_1
svmat ratios_sd_plus_1_plus_1, names(col)		

mat colnames ratios_sd_minus_1_minus_1 = SCC_sd_m_1_m_1 gamma_sd_m_1_m_1 cap_50HZ_sd_m_1_m_1 cap_50HZ_alloc_sd_m_1_m_1 cap_Amprion_sd_m_1_m_1 cap_Amprion_alloc_sd_m_1_m_1 cap_TenneT_sd_m_1_m_1 cap_TenneT_alloc_sd_m_1_m_1 cap_Transnet_sd_m_1_m_1 cap_Transnet_alloc_sd_m_1_m_1 v_bench_sd_m_1_m_1 v_alloc_sd_m_1_m_1 /// 
v_50HZ_sd_m_1_m_1 v_50HZ_alloc_sd_m_1_m_1 v_Amprion_sd_m_1_m_1 v_Amprion_alloc_sd_m_1_m_1 v_TenneT_sd_m_1_m_1 v_TenneT_alloc_sd_m_1_m_1 v_Transnet_sd_m_1_m_1 v_Transnet_alloc_sd_m_1_m_1
svmat ratios_sd_minus_1_minus_1, names(col)		

mat colnames ratios_sd_plus_1_minus_1 = SCC_sd_p_1_m_1 gamma_sd_p_1_m_1 cap_50HZ_sd_p_1_m_1 cap_50HZ_alloc_sd_p_1_m_1 cap_Amprion_sd_p_1_m_1 cap_Amprion_alloc_sd_p_1_m_1 cap_TenneT_sd_p_1_m_1 cap_TenneT_alloc_sd_p_1_m_1 cap_Transnet_sd_p_1_m_1 cap_Transnet_alloc_sd_p_1_m_1 v_bench_sd_p_1_m_1 v_alloc_sd_p_1_m_1 /// 
v_50HZ_sd_p_1_m_1 v_50HZ_alloc_sd_p_1_m_1 v_Amprion_sd_p_1_m_1 v_Amprion_alloc_sd_p_1_m_1 v_TenneT_sd_p_1_m_1 v_TenneT_alloc_sd_p_1_m_1 v_Transnet_sd_p_1_m_1 v_Transnet_alloc_sd_p_1_m_1
svmat ratios_sd_plus_1_minus_1, names(col)		

mat colnames ratios_sd_minus_1_plus_1 = SCC_sd_m_1_p_1 gamma_sd_m_1_p_1 cap_50HZ_sd_m_1_p_1 cap_50HZ_alloc_sd_m_1_p_1 cap_Amprion_sd_m_1_p_1 cap_Amprion_alloc_sd_m_1_p_1 cap_TenneT_sd_m_1_p_1 cap_TenneT_alloc_sd_m_1_p_1 cap_Transnet_sd_m_1_p_1 cap_Transnet_alloc_sd_m_1_p_1 v_bench_sd_m_1_p_1 v_alloc_sd_m_1_p_1 /// 
v_50HZ_sd_m_1_p_1 v_50HZ_alloc_sd_m_1_p_1 v_Amprion_sd_m_1_p_1 v_Amprion_alloc_sd_m_1_p_1 v_TenneT_sd_m_1_p_1 v_TenneT_alloc_sd_m_1_p_1 v_Transnet_sd_m_1_p_1 v_Transnet_alloc_sd_m_1_p_1
svmat ratios_sd_minus_1_plus_1, names(col)		


cap log close 
cap erase temp_v4.dta


// STORE MAIN VALUES
preserve
keep SCC-v_Transnet_alloc_sd_m_1_p_1 
drop if missing(SCC)
save "PVoutput_analysis3_10KW_v4_lpollution_SCC1.dta", replace
restore 

 

*******************************************************
// FIGURE D.4 (Panels b - d)

use "PVoutput_analysis3_10KW_v4_lpollution_SCC1.dta", clear


** Reallocation: SHARE of total capacity
gen cap_tot = cap_50HZ_alloc + cap_Amprion_alloc + cap_TenneT_alloc + cap_Transnet_alloc
gen alpha_50Hertz = cap_50HZ_alloc / cap_tot
gen alpha_Amprion = cap_Amprion_alloc  / cap_tot
gen alpha_TenneT = cap_TenneT_alloc  / cap_tot
gen alpha_TransnetBW = cap_Transnet_alloc / cap_tot

label var alpha_50Hertz "50Hertz"
label var alpha_Amprion "Amprion"
label var alpha_TenneT "TenneT"
label var alpha_TransnetBW "TransnetBW"
label var gamma "{&gamma}"


local solarcap_10kw_50Hertz   536.55963
local solarcap_10kw_Amprion   1779.9152
local solarcap_10kw_TenneT  2293.709
local solarcap_10kw_TransnetBW  1117.9264

local S 5728.1102

//qui sum gamma
//replace gamma = 0 if _n == `r(N)' + 1
cap drop dup
expand 2 if gamma == .5, gen(dup)
replace gamma = 0 if dup ==1

foreach v of varlist SCC cap_50HZ- alpha_TransnetBW {
	replace `v' = . if dup ==1
}
drop dup

sort gamma
gen beta_50Hertz = `solarcap_10kw_50Hertz' / `S' // if _n == `r(N)' + 1
gen beta_Amprion = `solarcap_10kw_Amprion' / `S' // if _n == `r(N)' + 1
gen beta_TenneT = `solarcap_10kw_TenneT' / `S' // if _n == `r(N)' + 1
gen beta_TransnetBW = `solarcap_10kw_TransnetBW' / `S' // if _n == `r(N)' + 1


label var beta_50Hertz "50Hertz"
label var beta_Amprion "Amprion"
label var beta_TenneT "TenneT"
label var beta_TransnetBW "TransnetBW"

sort gamma



tw scatter beta_50Hertz gamma if gamma == 0, ms(+) lw(medthick) col(ltblue) || ///
	scatter beta_Amprion gamma if gamma == 0, ms(d) lw(medthick) lp(-) col(cyan)  || ///
	scatter beta_TenneT gamma if gamma == 0, ms(o) lw(medthick) lp(.-) col(navy)  || ///
	scatter beta_TransnetBW gamma if gamma == 0, ms(t) lw(medthick) lp(_) col(teal) || ///
	line alpha_50Hertz gamma, lw(medthick) col(ltblue) || ///
	line alpha_Amprion gamma, lw(medthick) lp(-) col(cyan) || ///
	line alpha_TenneT gamma, lw(medthick) lp(.-) col(navy) || ///
	line alpha_TransnetBW gamma, lw(medthick) lp(_) col(teal) ylabel( , grid) ytitle("Ratio of solar capacity in TSO" "relative to total solar capacity") xtitle(`"{fontface "Times":{&gamma}}"') ///
  ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(alphas, replace)
graph export "../figures/ratios_alphas_lpollution_v4.pdf", replace 


*Total Gains as function of gamma
gen ratios = 100 * (v_alloc / v_bench - 1)

// THIS IS FOR THE BANDS
gen ratios_sd_plus_1_plus_1 = 100 * (v_alloc_sd_p_1_p_1 / v_bench_sd_p_1_p_1 - 1)
gen ratios_sd_minus_1_minus_1 = 100 * (v_alloc_sd_m_1_m_1 / v_bench_sd_m_1_m_1 - 1)
gen ratios_sd_minus_1_plus_1 = 100 * (v_alloc_sd_m_1_p_1 / v_bench_sd_m_1_p_1 - 1)
gen ratios_sd_plus_1_minus_1 = 100 * (v_alloc_sd_p_1_m_1 / v_bench_sd_p_1_m_1 - 1)
// MAKE SURE THAT COLUMN "MINUS_1" CONTAINS VALUES THAT ARE LOWER THAN THOSE IN COLUMN "PLUS_1":
egen ratios_sd_min = rowmin(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)
egen ratios_sd_max = rowmax(ratios_sd_plus_1_plus_1 ratios_sd_minus_1_minus_1 ratios_sd_minus_1_plus_1 ratios_sd_plus_1_minus_1)

gen ratios_AS = 100 * (v_AS_alloc / (v_AS_50HZ + v_AS_Amp + v_AS_TenneT + v_AS_Trans) - 1)		// PROPORTION OF ANC.SERV. VALUE AT THIS GAMA RELATIVE TO THE ANC.SERV. VALUE IN THE BASELINE
gen ratios_emissions = 100 * (v_ME_alloc / (v_ME_50HZ + v_ME_Amp + v_ME_TenneT + v_ME_Trans) - 1)
gen ratios_emissions_local = 100 * (v_ME_locpol_alloc / (v_ME_locpol_50HZ + v_ME_locpol_Amp + v_ME_locpol_TenneT + v_ME_locpol_Trans) - 1)
gen ratios_prodcosts = 100 * (v_MC_alloc / (v_MC_50HZ + v_MC_Amp + v_MC_TenneT + v_MC_Trans) - 1)


*Gains (by TSO) as function of gamma
tw 	line ratios_emissions gamma, lw(medthick) col(navy) || ///
	line ratios_emissions_local gamma, lw(medthick) lp(.-) col(dkgreen) || ///
	line ratios_prodcosts gamma, lw(medthick) lp(-) col(ltblue) || ///
	line ratios_AS  gamma, lw(medthick) lp(_) col(cyan) ///
	legend(lab(1 "Value of CO2 emissions") lab(2 "Value of co-pollutants")  lab(3 "Production costs") lab(4 "Ancillary service costs")) ylabel( , grid) ytitle("% change relative to baseline") xtitle(`"{fontface "Times":{&gamma}}"') ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(ratios_TSO_with_bands, replace)
graph export "../figures/ratios_decomposition_rel_baseline_lpollution_v4.pdf", replace 


*Decomposition of gains (total)
gen share_AS = 100 * v_AS_alloc / v_alloc
gen share_MC = 100 * v_MC_alloc / v_alloc
gen share_ME = 100 * v_ME_alloc / v_alloc
gen share_locpol_ME = 100 * v_ME_locpol_alloc / v_alloc


tw 	line share_ME gamma, lw(medthick) col(navy) || ///
	line share_locpol_ME gamma, lw(medthick) lp(.-) col(dkgreen) || ///
	line share_MC gamma, lw(medthick) lp(-) col(ltblue) || ///
	line share_AS gamma, lw(medthick) lp(_) col(cyan) ///
	legend(lab(1 "Value of CO2 emissions") lab(2 "Value of co-pollutants") lab(3 "Production costs") lab(4 "Ancillary service costs")) ylabel( , grid) ytitle("% of gains") xtitle(`"{fontface "Times":{&gamma}}"') ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(ratios_decomp, replace)
graph export "../figures/ratios_decomposition_lpollution_v4.pdf", replace 
*********************************************************************
